GDAS012-折衷t检验之limma分析微阵列数据

发表于 | 分类于
| 字数统计: 1,713 | 阅读时长 ≈ 8

简单t检验

在这一节中，我们将展示一下使用简单t检验和使用limma层次模型进行差异分析的区别。limma的相关参考文献已经列到参考资料中。

在这一节中，我们还民法了执行limma分析的基本步骤。需要注意的是，limma的功能非常强大，它有数百页的文档，我们不能详细地介绍这个包以及函数，有兴趣的同学可以直接查看它的文档。

spike-in数据集

我们先加载一批由Affymetrix生成的均一化(normalized)后的spike-in数据，如下所示：

1
2
3
4
# biocLite("SpikeInSubset")
library(SpikeInSubset)
data(rma95)
fac <- factor(rep(1:2,each=3))

简单来说，这里一共有16个mRNA的物种，它们都已经使用了固定的浓度，并混合起来使用了hgu95芯片进行了检验。这个子集中对于每个mRNA物种，ExpressionSet中的第一个三元组(trio)和第二个三元组有着固定的不同值。通过pData就能看到，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
pData(rma95)
##                  37777_at 684_at 1597_at 38734_at 39058_at 36311_at
## 1521a99hpp_av06      0.00   0.25     0.5        1        2        4
## 1532a99hpp_av04      0.00   0.25     0.5        1        2        4
## 2353a99hpp_av08      0.00   0.25     0.5        1        2        4
## 1521b99hpp_av06      0.25   0.50     1.0        2        4        8
## 1532b99hpp_av04      0.25   0.50     1.0        2        4        8
## 2353b99hpp_av08r     0.25   0.50     1.0        2        4        8
##                  36889_at 1024_at 36202_at 36085_at 40322_at 407_at
## 1521a99hpp_av06         8      16       32       64      128   0.00
## 1532a99hpp_av04         8      16       32       64      128   0.00
## 2353a99hpp_av08         8      16       32       64      128   0.00
## 1521b99hpp_av06        16      32       64      128      256   0.25
## 1532b99hpp_av04        16      32       64      128      256   0.25
## 2353b99hpp_av08r       16      32       64      128      256   0.25
##                  1091_at 1708_at 33818_at 546_at
## 1521a99hpp_av06      512    1024      256     32
## 1532a99hpp_av04      512    1024      256     32
## 2353a99hpp_av08      512    1024      256     32
## 1521b99hpp_av06     1024       0      512     64
## 1532b99hpp_av04     1024       0      512     64
## 2353b99hpp_av08r    1024       0      512     64

如果要了解这批数据这样设计的原理，现在我们来看一下下面的四张spiked mRNA，原文如下：

To get a feel for the response of the array quantifications to this design, consider the following plots for four of the spiked mRNAs.

1
2
3
4
5
6
7
8
par(mfrow=c(2,2))
for (i in 1:4) {
  spg = names(pData(rma95))
  plot(1:6, exprs(rma95)[spg[i+6],], main=spg[i+6], ylab="RMA",
    xlab="nominal", axes=FALSE)
  axis(2)
  axis(1, at=1:6, labels=pData(rma95)[[spg[i+6]]])
}

plot of chunk lkpl

由于RMA(robust multi-array average,鲁棒多阵列平均 )是以log2转换为定量标准的，因此在归一化的检测值中观察到的差异是比较合理的。但是，在每个设计好的浓度内部，存在着相当大的变异。

rowttests

我们现在使用 genefilter 包中的 rowttests 函数来做一个简单的t检验，如下所示： 

1
2
library(genefilter)
rtt <- rowttests(exprs(rma95),fac)

我们将根据p值，均值，以及特征值是否是spike-in值来设置颜色，如下所示：

1
2
3
mask <- with(rtt, abs(dm) < .2 & p.value < .01)
spike <- rownames(rma95) %in% colnames(pData(rma95))
cols <- ifelse(mask,"red",ifelse(spike,"dodgerblue","black"))

我们现在使用上面定义的颜色来绘制结果。我们将dm乘以-1，因为国我们感兴趣的是，第二组与第一组的差异（这是由lm与limma包默认情况下使用的差异）。spike-in基因是蓝色的，它们大多数有着较小的p值与较大的均值差异。红点表示较小p值的基因，同时也有着较小的均值差异。我们将会看到这些点使用limma后的变化，如下所示：

1
2
3
4
5
with(rtt, plot(-dm, -log10(p.value), cex=.8, pch=16,
     xlim=c(-1,1), ylim=c(0,5),
     xlab="difference in means",
     col=cols))
abline(h=2,v=c(-.2,.2), lty=2)

Note that the red genes have mostly low estimates of standard deviation.

1
2
3
4
rtt$s <- apply(exprs(rma95), 1, function(row) sqrt(.5 * (var(row[1:3]) + var(row[4:6]))))
with(rtt, plot(s, -log10(p.value), cex=.8, pch=16,
              log="x",xlab="estimate of standard deviation",
              col=cols))

plot of chunk unnamed-chunk-5

limma步骤

以下是执行基本limma分析的三个步骤。我们指定coef=2，因为我感兴趣的是组与组之间的差异，而不是截距(intercept)，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
library(limma)
options(digits=3)
fit <- lmFit(rma95, design=model.matrix(~ fac))  # step 1 least squares estimates
colnames(coef(fit))     
## [1] "(Intercept)" "fac2"
fit <- eBayes(fit)           # step 2 moderate the t statistics
tt <- topTable(fit, coef=2)  # step 3 report
tt
##           logFC AveExpr      t  P.Value adj.P.Val    B
## 1708_at  -7.061    7.95 -73.53 7.82e-17  9.87e-13 8.65
## 36202_at  0.853    9.37   9.98 4.94e-07  3.12e-03 4.59
## 36311_at  0.832    8.56   8.36 3.02e-06  1.27e-02 3.57
## 33264_at  0.712    4.92   7.43 9.67e-06  2.71e-02 2.84
## 32660_at  0.655    8.68   7.36 1.07e-05  2.71e-02 2.77
## 38734_at  0.747    6.26   7.19 1.35e-05  2.83e-02 2.62
## 1024_at   0.843    9.70   6.73 2.50e-05  4.40e-02 2.20
## 36085_at  0.645   12.19   6.65 2.79e-05  4.40e-02 2.12
## 33818_at  0.532   12.29   6.45 3.70e-05  5.19e-02 1.92
## 39058_at  0.609    7.53   6.28 4.77e-05  5.69e-02 1.74

topTable 会返回你定义的值排序的前几个基因。你还可以使用topTable通过指定none来要求返回所有的值。需要注意的是，有一列会给出校正后的每个基因的p值。默认情况下，通过p.adjust函数中的Benjamini-Hochberg方法用于校正p值，如下所示：

1
2
3
4
# ?topTable
dim(topTable(fit, coef=2, number=Inf, sort.by="none"))
## [1] 12626     6
# ?p.adjust

在这里，我们比较一下使用limma分析的结果与以前得到的火山图结果。注意，红点现在都位于-log10(p.value)=2之下。此外，表示实际差异的蓝点的p值比以前更高，如下所示：

1
2
3
4
5
limmares <- data.frame(dm=coef(fit)[,"fac2"], p.value=fit$p.value[,"fac2"])
with(limmares, plot(dm, -log10(p.value),cex=.8, pch=16,
     col=cols,xlab="difference in means",
     xlim=c(-1,1), ylim=c(0,5)))
abline(h=2,v=c(-.2,.2), lty=2)

最后我们来构建一张图，用于展示limma是如何将方差估计缩小到一个共同值(common value)的，这能消除由于方差估计过低而导致的假阳性。

在这里，我们为不同的方差估计值设计40个小区间(bin)，每个区间中落入一个基因。我们会删除那些不含有任何基因的区间，原文如下：

Here we pick, for each of 40 bins of different variance estimates, a single gene which falls in that bin. We remove bins which do not have any such genes.

1
2
3
4
n <- 40
qs <- seq(from=0,to=.2,length=n)
idx <- sapply(seq_len(n),function(i) which(as.integer(cut(rtt$s^2,qs)) == i)[1])
idx <- idx[!is.na(idx)]

现在绘制一条直线，这个直线展示了那些基因最初的估计方差到运行了limma之后的估计方差，如下所示：

1
2
3
4
5
6
par(mar=c(5,5,2,2))
plot(1,1,xlim=c(0,.21),ylim=c(0,1),type="n",
     xlab="variance estimates",ylab="",yaxt="n")
axis(2,at=c(.1,.9),c("before","after"),las=2)
segments((rtt$s^2)[idx],rep(.1,n),
         fit$s2.post[idx],rep(.9,n))

plot of chunk unnamed-chunk-10

关于limma校正的一些层级模型资料，可以参考以前的笔记《DALS019-统计模型2-贝叶斯分布与层次分析》。

参考资料

  1. Using limma for microarray analysis
  2. Smyth GK, “Linear models and empirical bayes methods for assessing differential expression in microarray experiments”. Stat Appl Genet Mol Biol. 2004 http://www.ncbi.nlm.nih.gov/pubmed/16646809